rm(list = ls())
library(tidyr)
library(readr)
library(dplyr)
library(readxl)
library(tidyr)
library(arm)
library(lme4)
library(MASS)
library(texreg)
load("Data/paperdata.RData") #1993-2006

source("Code/coefplot.R")
summary(paperdata)

names(paperdata)
# direct: x/p
# indirect: u, o, l,i, f,t,m,w,y

paperdata <- paperdata %>% 
  dplyr::mutate(ext_direct = ifelse((ext_x == 0 & ext_p == 0), 0, 1))

paperdata <- paperdata %>% 
  dplyr::mutate(ext_indirect = ifelse((ext_y == 0 & ext_w == 0 & ext_m == 0 & 
                                         ext_t == 0 & ext_f == 0 & ext_i == 0 &
                                         ext_l == 0 & ext_o == 0 & ext_u == 0), 0, 1))

## indrect-gov
paperdata <- paperdata %>% 
  dplyr::mutate(ext_gov_indirect = ifelse((ext_y_gov == 0 & ext_w_gov == 0 & ext_m_gov == 0 & 
                                             ext_t_gov == 0 & ext_f_gov == 0 & ext_i_gov == 0 &
                                             ext_l_gov == 0 & ext_o_gov == 0 & ext_u_gov == 0), 0, 1))
## dirct-rebel
paperdata <- paperdata %>% 
  dplyr::mutate(ext_rebel_indirect = ifelse((ext_y_rebel == 0 & ext_w_rebel == 0 & ext_m_rebel == 0 & 
                                               ext_t_rebel == 0 & ext_f_rebel == 0 & ext_i_rebel == 0 &
                                               ext_l_rebel == 0 & ext_o_rebel == 0 & ext_u_rebel == 0), 0, 1))

paperdata <- paperdata %>% dplyr::select(ext_direct,ext_indirect, ext_x, ext_p, everything())

## DV: 
# ext_count
#ext_x_count
#ext_p
# ext_y
#ext_y_coun
#####
# library(spduration)
# 
# #transform to spdur format
# paperdata <- as.data.frame(paperdata)
# paperdata <- add_duration(paperdata, "ext_direct", unitID = "dyad_id", tID ="year", freq = "year", ongoing = FALSE)
# 
# ## create no-alliance years polynomial terms
# paperdata<- paperdata  %>% 
#   dplyr::mutate(t_yrdirect = duration,
#                 t_yrsqdirect = duration^2,
#                 t_yrcubdirect = duration^3)


dv <- "ext_direct"

## controls
ctrs <- c("polity2_lag1","gdp_log_lag1", "military_log_lag1", "population_log_lag1",
          "bd_best_lag1", "deaths_civilians_lag1", "factor(conflict_id)","year")
ivs_1 <- c("troops_backward_lag1")
ivs_2 <- c("police_backward_lag1")
ivs_3 <- c("observers_backward_lag1")
ivs_4 <- c("totals_backward_lag1")


f_1 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_1, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))
f_2 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_2, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))
f_3 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_3, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))
f_4 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_4, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))
ivs <- c("troops_backward_lag1","police_backward_lag1","observers_backward_lag1","totals_backward_lag1")

library(arm)
## modeling
model_ext_sup_direct <- list()
for (i in 1:length(ivs)){
  model_ext_sup_direct[[i]] <- glm(eval(parse(text = paste0('f_', i))), 
                                   data = paperdata, family = binomial(link = "logit"))
}

# for (i in 1:length(ivs)){
#   model_ext_sup_direct[[i]] <- bayesglm(eval(parse(text = paste0('f_', i))), 
#                                    data = paperdata, family = binomial(link = "logit"),prior.scale=2.5, prior.df=1)
# }


dv <- "ext_indirect"

## controls
ctrs <- c("polity2_lag1","gdp_log_lag1", "military_log_lag1", "population_log_lag1",
          "bd_best_lag1", "deaths_civilians_lag1", "factor(gwno_loc)", "factor(year)")
ivs_1 <- c("troops_backward_lag1")
ivs_2 <- c("police_backward_lag1")
ivs_3 <- c("observers_backward_lag1")
ivs_4 <- c("totals_backward_lag1")


f_1 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_1, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))
f_2 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_2, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))
f_3 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_3, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))
f_4 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_4, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))
ivs <- c("troops_backward_lag1","police_backward_lag1","observers_backward_lag1","totals_backward_lag1")


## modeling
model_ext_sup_indirect <- list()
for (i in 1:length(ivs)){
  model_ext_sup_indirect[[i]] <- glm(eval(parse(text = paste0('f_', i))), 
                                     data = paperdata, family = binomial(link = "logit"))
}


p1_1 <- marginaleffect_logit2(ModelResults =  model_ext_sup_direct  , n.sim = 1000, 
                              varname = c("troops_backward_lag1","police_backward_lag1",
                                          "observers_backward_lag1","totals_backward_lag1"),
                              data = paperdata,
                              dv = "ext_direct",
                              val1 = 0, val2 =1, clusterid = "gwno_loc") 


p1_2 <- marginaleffect_logit2(ModelResults = model_ext_sup_indirect, n.sim = 1000, 
                              varname = c("troops_backward_lag1","police_backward_lag1",
                                          "observers_backward_lag1","totals_backward_lag1"),
                              data = paperdata,
                              dv = "ext_indirect",
                              val1 = 0, val2 =1, clusterid = "gwno_loc") 




df_ext_sup_nonUN <- bind_rows(p1_1, p1_2)



p = ggplot(df_ext_sup_nonUN, aes(colour = dv)) + 
  geom_hline(yintercept = 0, lty = 2) +
  geom_linerange(aes(x = id, ymin = low,
                     ymax = high),
                 lwd = 1.5, position = position_dodge(width = 1/2)) +
  geom_pointrange(aes(x = id, y =median, ymin =low,
                      ymax =high, shape = dv),
                  fatten = 7,
                  lwd = 1.5, position = position_dodge(width = 1/2),
                  fill="white")+
  geom_text(aes(x = id +0.2 ,
                y = median,
                label = format(round(median, 2))),position = position_dodge(width = 1/2)) +
  theme_bw() + xlab("") + ylab("预测概率变化（边际效应）") + 
  theme(legend.position="bottom",
        legend.title=element_blank(),
        legend.text = element_text(colour="black", size=16),
        axis.text= element_text(size=16),
        text = element_text(size=16, family = 'KaiTi',face = "bold"),
        plot.title = element_text(hjust = .5, size = 16))+
  ggtitle("因变量：是否存在外部支持")+
  labs(tag = "图5a ")+
  theme(plot.tag.position = "bottom")+
  scale_colour_manual(values = c("black", "black"), labels = c("直接支持","间接支持"))+
  scale_shape_manual(values = c(15,23),labels = c("直接支持","间接支持"))+
  guides(color = guide_legend(nrow = 1))+
  scale_x_continuous(breaks = 1:4,
                     labels = parse(text =c( 
                       "维和部队",
                       "维和警察",
                       "维和观察员",
                       "所有类型")))
ggsave("Results/figures/fd_ext_sup_directindirect.pdf", width = 6, height = 4)
ggsave("Replication/figures/figure5-a.bmp", width = 6, height = 4)
ggsave("Replication/figures/figure5-a.pdf", width = 6, height = 4)



save(model_ext_sup_direct,model_ext_sup_indirect, file = "Results/directindirect.RData")



modSD = cluster_se(ModelResults = c(model_ext_sup_direct,model_ext_sup_indirect), data = paperdata, 
                   clusterid = "gwno_loc") 

modPvalue = cluster_pvalue(ModelResults = c(model_ext_sup_direct,model_ext_sup_indirect), data = paperdata, 
                           clusterid = "gwno_loc") 


screenreg(c(model_ext_sup_direct,model_ext_sup_indirect ),override.se =  modSD)
## table-1
htmlreg(c(model_ext_sup_direct,model_ext_sup_indirect ), file = "Replication/tables/directindirect.doc",
        stars = c(0.01, 0.05, 0.1),
        override.se =  modSD,
        override.pvalues = modPvalue,
        caption = "Logit回归结果：维和行动对外部支持直接与间接支持方式的影响",
        caption.above = TRUE,
        label = "tab:tab5",
        scalebox = 0.7,
        custom.header = list("直接支持" = 1:4,
                             "间接支持" = 5:8),
        # custom.model.names = c("M1:troop", "M2:police", "M3:observers", "M4:any",
        #                        "M5:troop", "M6:police", "M7:observers", "M8:any",
        #                        "M9:troop", "M10:police", "M11:observers", "M12:any"),
        custom.coef.map = list("troops_backward_lag1" = "维和部队",
                               "police_backward_lag1" = "维和警察",
                               "observers_backward_lag1" ="维和观察员",
                               "totals_backward_lag1" = "所有类型",
                               "polity2_lag1" = "政体类型",
                               "gdp_log_lag1" = "国民生产总值（对数）",
                               "military_log_lag1" = "军队人数（对数）",
                               "population_log_lag1" = "人口总数(对数)",
                               "bd_best_lag1"= "战斗相关死亡人数",
                               "deaths_civilians_lag1" = "平民死亡人数",
                               "year" = "年份"),
        custom.gof.rows = list("国家固定效应" = c("否", "否", "否", "否",
                                            "是", "是", "是", "是"),
                               "冲突固定效应" = c("是", "是", "是", "是",
                                            "否", "否", "否", "否"),
                               "年份固定效应" = c("否", "否", "否", "否",
                                            "是", "是", "是", "是")),
        digits = 3)














### 间接支持-un

dv <- "ext_gov_indirect"

## controls
ctrs <- c("polity2_lag1","gdp_log_lag1", "military_log_lag1", "population_log_lag1",
          "bd_best_lag1", "deaths_civilians_lag1", "factor(gwno_loc)", "factor(year)")
ivs_1 <- c("troops_backward_un_lag1")
ivs_2 <- c("police_backward_un_lag1")
ivs_3 <- c("observers_backward_un_lag1")
ivs_4 <- c("totals_backward_un_lag1")


f_1 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_1, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))
f_2 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_2, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))
f_3 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_3, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))
f_4 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_4, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))
ivs <- c("troops_backward_un_lag1","police_backward_un_lag1","observers_backward_un_lag1","totals_backward_un_lag1")


## modeling
model_ext_gov_indirect_un <- list()
for (i in 1:length(ivs)){
  model_ext_gov_indirect_un[[i]] <- glm(eval(parse(text = paste0('f_', i))), 
                                        data = paperdata, family = binomial(link = "logit"))
}

dv <- "ext_rebel_indirect"

## controls
ctrs <- c("polity2_lag1","gdp_log_lag1", "military_log_lag1", "population_log_lag1",
          "bd_best_lag1", "deaths_civilians_lag1", "factor(gwno_loc)", "factor(year)")
ivs_1 <- c("troops_backward_un_lag1")
ivs_2 <- c("police_backward_un_lag1")
ivs_3 <- c("observers_backward_un_lag1")
ivs_4 <- c("totals_backward_un_lag1")


f_1 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_1, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))
f_2 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_2, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))
f_3 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_3, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))
f_4 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_4, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))
ivs <- c("troops_backward_un_lag1","police_backward_un_lag1","observers_backward_un_lag1","totals_backward_un_lag1")


## modeling
model_ext_rebel_indirect_un <- list()
for (i in 1:length(ivs)){
  model_ext_rebel_indirect_un[[i]] <- glm(eval(parse(text = paste0('f_', i))), 
                                          data = paperdata, family = binomial(link = "logit"))
}




#
p1_1 <- marginaleffect_logit2(ModelResults =  model_ext_gov_indirect_un   , n.sim = 1000, 
                              varname = c("troops_backward_un_lag1","police_backward_un_lag1",
                                          "observers_backward_un_lag1","totals_backward_un_lag1"),
                              data = paperdata,
                              dv = "ext_gov_indirect",
                              val1 = 0, val2 =1, clusterid = "gwno_loc") 


p1_2 <- marginaleffect_logit2(ModelResults = model_ext_rebel_indirect_un, n.sim = 1000, 
                              varname = c("troops_backward_un_lag1","police_backward_un_lag1",
                                          "observers_backward_un_lag1","totals_backward_un_lag1"),
                              data = paperdata,
                              dv = "ext_rebel_indirect",
                              val1 = 0, val2 =1, clusterid = "gwno_loc") 




df_ext_indrect_un <- bind_rows(p1_1, p1_2)



p = ggplot(df_ext_indrect_un, aes(colour = dv)) + 
  geom_hline(yintercept = 0, lty = 2) +
  geom_linerange(aes(x = id, ymin = low,
                     ymax = high),
                 lwd = 1.5, position = position_dodge(width = 1/2)) +
  geom_pointrange(aes(x = id, y =median, ymin =low,
                      ymax =high, shape = dv),
                  fatten = 7,
                  lwd = 1.5, position = position_dodge(width = 1/2),
                  fill="white")+
  geom_text(aes(x = id +0.2 ,
                y = median,
                label = format(round(median, 2))),position = position_dodge(width = 1/2)) +
  theme_bw() + xlab("") + ylab("预测概率变化（边际效应）") + 
  theme(legend.position="bottom",
        legend.title=element_blank(),
        legend.text = element_text(colour="black", size=16),
        axis.text= element_text(size=16),
        text = element_text(size=16, family = 'KaiTi',face = "bold"),
        plot.title = element_text(hjust = .5, size = 16))+
  ggtitle("因变量：是否存在外部间接支持(联合国维和)")+
  labs(tag = "图5a ")+
  theme(plot.tag.position = "bottom")+
  scale_colour_manual(values = c("black", "black"), labels = c("支持政府","支持叛乱组织"))+
  scale_shape_manual(values = c(15,23),labels = c("支持政府","支持叛乱组织"))+
  guides(color = guide_legend(nrow = 1))+
  scale_x_continuous(breaks = 1:4,
                     labels = parse(text =c( 
                       "维和部队",
                       "维和警察",
                       "维和观察员",
                       "所有类型")))
ggsave("Results/figures/fd_ext_un_indirect.pdf", width = 6, height = 4)
ggsave("Replication/figures/figure5-b.bmp", width = 6, height = 4)
ggsave("Replication/figures/figure5-b.pdf", width = 6, height = 4)





save(model_ext_gov_indirect_un,model_ext_rebel_indirect_un, file = "Results/ext_sup_UN_indirect.RData")
modSD = cluster_se(ModelResults = c(model_ext_gov_indirect_un,model_ext_rebel_indirect_un), data = paperdata, 
                   clusterid = "gwno_loc") 
modPvalue = cluster_pvalue(ModelResults = c(model_ext_gov_indirect_un,model_ext_rebel_indirect_un), data = paperdata, 
                           clusterid = "gwno_loc") 


screenreg(c(model_ext_gov_indirect_un,model_ext_rebel_indirect_un))

htmlreg(c(model_ext_gov_indirect_un,model_ext_rebel_indirect_un), file = "Results/tables/ext_sup_UN_indirect.doc",
        stars = c(0.01, 0.05, 0.1),
        override.se =  modSD,
        override.pvalues = modPvalue,
        caption = "Logit回归结果：联合国维和行动对外部间接支持的影响",
        caption.above = TRUE,
        label = "tab:tab5",
        scalebox = 0.7,
        custom.header = list("外部政府" = 1:4,
                             "支持反叛组织" = 5:8),
        # custom.model.names = c("M1:troop", "M2:police", "M3:observers", "M4:any",
        #                        "M5:troop", "M6:police", "M7:observers", "M8:any",
        #                        "M9:troop", "M10:police", "M11:observers", "M12:any"),
        custom.coef.map = list("troops_backward_un_lag1" = "维和部队",
                               "police_backward_un_lag1" = "维和警察",
                               "observers_backward_un_lag1" ="维和观察员",
                               "totals_backward_un_lag1" = "所有类型",
                               "polity2_lag1" = "政体类型",
                               "gdp_log_lag1" = "国民生产总值（对数）",
                               "military_log_lag1" = "军队人数（对数）",
                               "population_log_lag1" = "人口总数(对数)",
                               "bd_best_lag1"= "战斗相关死亡人数",
                               "deaths_civilians_lag1" = "平民死亡人数"),
        custom.gof.rows = list("国家固定效应" = c("是", "是", "是", "是",
                                            "是", "是", "是", "是"),
                               "年份固定效应" = c("是", "是", "是", "是",
                                            "是", "是", "是", "是")),
        digits = 3)


#### non-Un
dv <- "ext_gov_indirect"

## controls
ctrs <- c("polity2_lag1","gdp_log_lag1", "military_log_lag1", "population_log_lag1",
          "bd_best_lag1", "deaths_civilians_lag1", "factor(gwno_loc)", "factor(year)")
ivs_1 <- c("troops_backward_non_lag1")
ivs_2 <- c("police_backward_non_lag1")
ivs_3 <- c("observers_backward_non_lag1")
ivs_4 <- c("totals_backward_non_lag1")


f_1 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_1, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))
f_2 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_2, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))
f_3 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_3, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))
f_4 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_4, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))
ivs <- c("troops_backward_non_lag1","police_backward_non_lag1",
         "observers_backward_non_lag1","totals_backward_non_lag1")


## modeling
model_ext_gov_indirect_non <- list()
for (i in 1:length(ivs)){
  model_ext_gov_indirect_non[[i]] <- glm(eval(parse(text = paste0('f_', i))), 
                                         data = paperdata, family = binomial(link = "logit"))
}

dv <- "ext_rebel_indirect"

## controls
ctrs <- c("polity2_lag1","gdp_log_lag1", "military_log_lag1", "population_log_lag1",
          "bd_best_lag1", "deaths_civilians_lag1", "factor(gwno_loc)", "factor(year)")
ivs_1 <- c("troops_backward_non_lag1")
ivs_2 <- c("police_backward_non_lag1")
ivs_3 <- c("observers_backward_non_lag1")
ivs_4 <- c("totals_backward_non_lag1")


f_1 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_1, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))
f_2 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_2, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))
f_3 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_3, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))
f_4 <- formula(paste0(dv, ' ~ ', paste(paste(ivs_4, collapse=' + '), paste(ctrs, collapse=' + '), sep =' + ')))
ivs <- c("troops_backward_non_lag1","police_backward_non_lag1",
         "observers_backward_non_lag1","totals_backward_non_lag1")


## modeling
model_ext_rebel_indirect_non <- list()
for (i in 1:length(ivs)){
  model_ext_rebel_indirect_non[[i]] <- glm(eval(parse(text = paste0('f_', i))), 
                                           data = paperdata, family = binomial(link = "logit"))
}




#
p1_1 <- marginaleffect_logit2(ModelResults =  model_ext_gov_indirect_non   , n.sim = 1000, 
                              varname = c("troops_backward_non_lag1","police_backward_non_lag1",
                                          "observers_backward_non_lag1","totals_backward_non_lag1"),
                              data = paperdata,
                              dv = "ext_gov_indirect",
                              val1 = 0, val2 =1, clusterid = "gwno_loc") 


p1_2 <- marginaleffect_logit2(ModelResults = model_ext_rebel_indirect_non, n.sim = 1000, 
                              varname = c("troops_backward_non_lag1","police_backward_non_lag1",
                                          "observers_backward_non_lag1","totals_backward_non_lag1"),
                              data = paperdata,
                              dv = "ext_rebel_indirect",
                              val1 = 0, val2 =1, clusterid = "gwno_loc") 




df_ext_indrect_non <- bind_rows(p1_1, p1_2)



p = ggplot(df_ext_indrect_non, aes(colour = dv)) + 
  geom_hline(yintercept = 0, lty = 2) +
  geom_linerange(aes(x = id, ymin = low,
                     ymax = high),
                 lwd = 1.5, position = position_dodge(width = 1/2)) +
  geom_pointrange(aes(x = id, y =median, ymin =low,
                      ymax =high, shape = dv),
                  fatten = 7,
                  lwd = 1.5, position = position_dodge(width = 1/2),
                  fill="white")+
  geom_text(aes(x = id +0.2 ,
                y = median,
                label = format(round(median, 2))),position = position_dodge(width = 1/2)) +
  theme_bw() + xlab("") + ylab("预测概率变化（边际效应）") + 
  theme(legend.position="bottom",
        legend.title=element_blank(),
        legend.text = element_text(colour="black", size=16),
        axis.text= element_text(size=16),
        text = element_text(size=16, family = 'KaiTi',face = "bold"),
        plot.title = element_text(hjust = .5, size = 16))+
  ggtitle("因变量：是否存在外部间接支持(非联合国维和)")+
  labs(tag = "图5b ")+
  theme(plot.tag.position = "bottom")+
  scale_colour_manual(values = c("black", "black"), labels = c("支持政府","支持叛乱组织"))+
  scale_shape_manual(values = c(15,23), labels = c("支持政府","支持叛乱组织"))+
  guides(color = guide_legend(nrow = 1))+
  scale_x_continuous(breaks = 1:4,
                     labels = parse(text =c( 
                       "维和部队",
                       "维和警察",
                       "维和观察员",
                       "所有类型")))
ggsave("Results/figures/fd_ext_non_indirect.pdf", width = 6, height = 4)
ggsave("Replication/figures/figure5-c.bmp", width = 6, height = 4)
ggsave("Replication/figures/figure5-c.pdf", width = 6, height = 4)


save(model_ext_gov_indirect_non,model_ext_rebel_indirect_non, file = "Results/ext_sup_nonUN_indirect.RData")
modSD = cluster_se(ModelResults = c(model_ext_gov_indirect_non,model_ext_rebel_indirect_non), data = paperdata, 
                   clusterid = "gwno_loc") 

modPvalue = cluster_pvalue(ModelResults = c(model_ext_gov_indirect_non,model_ext_rebel_indirect_non), data = paperdata, 
                           clusterid = "gwno_loc") 

screenreg(c(model_ext_gov_indirect_non,model_ext_rebel_indirect_non))

htmlreg(c(model_ext_gov_indirect_non,model_ext_rebel_indirect_non ), file = "Replication/tables/ext_sup_nonUN_indirect.doc",
        stars = c(0.01, 0.05, 0.1),
        override.se =  modSD,
        override.pvalues = modPvalue,
        caption = "Logit回归结果：非联合国组织维和行动对外部间接支持的影响",
        caption.above = TRUE,
        label = "tab:tab5",
        scalebox = 0.7,
        custom.header = list("支持政府" = 1:4,
                             "支持反叛组织" = 5:8),
        # custom.model.names = c("M1:troop", "M2:police", "M3:observers", "M4:any",
        #                        "M5:troop", "M6:police", "M7:observers", "M8:any",
        #                        "M9:troop", "M10:police", "M11:observers", "M12:any"),
        custom.coef.map = list("troops_backward_non_lag1" = "维和部队",
                               "police_backward_non_lag1" = "维和警察",
                               "observers_backward_non_lag1" ="维和观察员",
                               "totals_backward_non_lag1" = "所有类型",
                               "polity2_lag1" = "政体类型",
                               "gdp_log_lag1" = "国民生产总值（对数）",
                               "military_log_lag1" = "军队人数（对数）",
                               "population_log_lag1" = "人口总数(对数)",
                               "bd_best_lag1"= "战斗相关死亡人数",
                               "deaths_civilians_lag1" = "平民死亡人数"),
        custom.gof.rows = list("国家固定效应" = c("是", "是", "是", "是",
                                            "是", "是", "是", "是"),
                               "年份固定效应" = c("是", "是", "是", "是",
                                            "是", "是", "是", "是")),
        digits = 3)




#write_dta(paperdata, "Data/paperdata.dta")
